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ABSTRACT 

Cosmological shock waves result from supersonic flow motions induced by hi- 
erarchical clustering of nonlinear structures in the universe. These shocks govern 
the nature of cosmic plasma through thermalization of gas and acceleration of 
nonthermal, cosmic-ray (CR) particles. We study the statistics and energetics of 
shocks formed in cosmological simulations of a concordance ACDM universe, with 
a special emphasis on the effects of non-gravitational processes such as radiative 
cooling, photoionization/heating, and galactic superwind feedbacks. Adopting 
an improved model for gas thermalization and CR acceleration efficiencies based 
on nonlinear diffusive shock acceleration calculations, we then estimate the gas 
thermal energy and the CR energy dissipated at shocks through the history of 
the universe. Since shocks can serve as sites for generation of vorticity, we also 
examine the vorticity that should have been generated mostly at curved shocks in 
cosmological simulations. We find that the dynamics and energetics of shocks are 
governed primarily by the gravity of matter, so other non-gravitational processes 
do not affect significantly the global energy dissipation and vorticity generation 
at cosmological shocks. Our results reinforce scenarios in which the intraclus- 
ter medium and warm-hot intergalactic medium contain energetically significant 
populations of nonthermal particles and turbulent fiow motions. 

Subject headings: cosmic rays - large-scale structure of universe - methods: nu- 
merical - shock waves - turbulence 
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Introduction 



Astrophysical plasmas consist of botli thermal particles and nonthermal, cosmic-ray 
(CR) particles that are closely coupled with permeating magnetic fields and underlying 
turbulent fiows. In the interstellar medium (ISM) of our Galaxy, for example, an approx- 
imate energy equipartition among different c omponents see ms to have been established, 
i.e., Etherm ~ ^CR ~ £^_B ~ £turb ~ 1 cV cm~^ (jLongaiii Il994l ) ■ Understanding the complex 
network of physical interactions among these components constitutes one of fundamental 
problems in astrophysics. 

There is substantial observational evidence for the presence of nonthermal particles and 
magnetic fields in the large scale structure of the universe. A fair fraction of X-ray clusters 
have been observed in diffuse radio synchrotron emission, indi cating the presence of GeV 
CR electrons and nG fields in the intracluster medium (ICM) (iGiovannini fc Ferettill2000l ). 
Observations in EUV and hard X-ray have shown that some clusters possess excess radia- 
tion compared to what is expected from the hot, thermal X-ray emitting ICM, most likely 
produced by the inve rse-Compton scattering of cosmic background radiation (CBR) pho - 
tons by CR electrons (IFusco-Femiano et a/.lll999l : iBowyer et a/.lll999l : iBerghofer et a/.ll2000l ). 
Assuming energy equipartition between CR electrons and magnetic fiel ds, ecRe ^ £r 



0.01— O.leV cm ~ 10 10 etherm can be inferred in typical radio halos (iGovoni fc Feretti 



20041 ). If some of those CR electrons have been energized at shocks and/or by turbulence, 
the same process should have produced a greater CR proton po pulation. Considerin g the 
ratio of proton to electron numbers, K ~ 100, for Galactic CRs (IBeck fc Krausll2005l ). one 
can expect ecRp ~ 0.01 — O.letherm in radio halos. However, CR protons in the ICM have 
yet to be confirmed by the observation of 7-ray photons pro duced by inelastic collisions 
between CR protons and thermal protons (IReimer et a/.ll2003l ). Magnetic fields have been 
also directly observed with Faraday rotation measure (RM). In clusters of galaxies strong 
fields of a few /iG strength extending from core t o 500 kpc or further were inferred from 
RM observations (jClarke et a/.ll200ll : IClarke II2004I ). An upper limit of < /iG was imposed 
on the magnetic field stren gth in filaments and sheets, based the observed limit of the RMs 
of quasars outside clusters (lKronberglll994 iRyu et a/.lll998l ). 



Studies on turbulence and turbulent magnetic fields in the large scale structure of the 
universe have been recently launched too. XMM-Newton X-ray observations of the Coma 
cluster, which seems to be in a post-merger stag e, were analyzed in details to extract clues 
on turbulence in the ICM (jSchuecker et al\\2004 ). By analyzing pressure fiuctuations, it was 
shown that the turbulence is likely subsonic and consistent with Kolmogoroff turbulence. 
RM maps of clusters h ave been analyzed to find the power spe ctrum o f turbulent magnetic 
fields in a few clusters jMurgia et alhooi IVogt fc EnBlinlboOsh . While iMurgia et all tooj ) 
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reported a spectrum shallower than the Kolmogoroff spectrum, IVogt fc EnBliru (120051 ) argued 
that the spectrum could be consistent with the Kolmogoroff spectrum if it is bended at a 
few kpc scale. These studies suggest that as in the ISM, turbulence does exist in the ICM 
and may constitute an energetically non-negligible component. 

In galaxy cluster environments there are s everal possible sources of CRs, mag netic 
fields, and turbulence: jets from active galaxies (IKronberg et a/.ll2004l: iLi et alluOQw . ter- 



mination shock s of galactic winds driven by supernova explosions (Vok fc AtoyanI Il999l ) 



merger shocks (ISarazin 



19991 



shocks ( Loeb &: Waxmann 200ol: 



Gabici fc Blasi 



2OO3I: iFujita et a/.ll2003l ). structure formation 



Miniati et al. 



2001a 



b|), and motions of subcluster clumps 



and galaxies (ISubramanian et al\ l2006l ) . All of them have a potential to inject a similar 
amount of energies, i.e., E ^ 10^^ — 10^^ ergs into the ICM. Here we focus on shock scenar- 
ios. 

Astrophysical shocks are collisionless shocks that form in tenuous cosmic plasmas via col- 
lective electromagnetic interactions between gas particles and magnetic fields. They play key 
roles in governing the nature of cosmic plasmas: i.e., 1) shocks convert a part of the kinetic 
energy of bulk flow motions into thermal energy , 2) shocks accelerate CRs by diffusive shock 



accel eration (DSA) (IBlandford fc Ostrikerlll978l : lBlandford &: Eichler 



20011), and amplify magnetic fields by streaming CRs flBelll Il978l : iLucek fc Belli I20001). 3) 



1987 



Malkov &: Drurv 



1950 



shocks generate m agnetic fields via the Bierr aann battery mechanism (|Biermann 
Kulsrud a/1ll997h and the Weibel instability JWeibellE959 : iMedvedev et a/.ll2006h. and 41 
curve d shocks generate vorticity and ensuing turbulent flows (lBinneylll974l : lDavies &: Widrow 
2000h . 



In iRyu et al\ (120031 ) (Paper I), the properties of cosmological shock waves in the inter- 
galactic medium (IGM) and the energy dissipations into thermal and nonthermal components 
at those shocks were studied in a high-resolution, adiabatic (non-radiative), hydrodynamic 
simulation of a ACDM universe. They found that internal shocks with low Mach numbers 
of M < 4, which formed in the hot, previously shocked gas inside nonlinear structures, 
are responsible for most of the shock energy dissipation. Adopting a nonlinear DSA model 
for CR protons, it was shown that about 1/2 of the gas thermal energy dissipated at cos- 
molog ical shocks through the history of the universe could be stored as CRs. In a recent 
study, iPfrommer et al\ (120061 ) identified shocks and analyzed the statistics in smoothed par- 
ticle hydrodynamic (SPH) simulations of a ACDM universe, and found that their results 
are in good agreement with those of Paper I. While internal shocks with lower Mach num- 
bers are energetically dominant, external accretions shocks with hi gher Mach numbers can 
serve as possible acceleration sites f or high energy cosmic rays (IKang et al\ Il996l . Il997l : 
Ostrowski &: Siemieniec-Oziebld 120021 ). It was shown that CR ions could be accelerated up 
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to ~ Z X 10 at cosmological shocks, where Z is the charge of ions (llnoue et a/.l 120071 ) . 



Ryu et al\ (120071 ) (Paper II) analyzed the distribution of vorticity, which should have 



been generated mostly at cosmological shock waves, in the same simulation of a ACDM 
universe as in Paper I, and studied its implication on turbulence and turbulence dynamo. 
Inside nonlinear structures, vorticity was found to be large enough that the turn-over time, 
w hich is defined as the inverse of vorticity, is shorter than the age of the universe. Based on 



it IRyu et a/.l (120071 ) argued that turbulence should have been developed in those structures 
and estimated the strength of the magnetic field grown by the turbulence. 

In this paper, we study cosmological shock waves in a new set of hydrodynamic sim- 
ulations of large structure formation in a concordance ACDM universe: an adiabatic (non- 
radiative) simulation which is similar to that considered in Paper I, and two additional 
simulations which include various non-gravitational processes (see the next section for de- 
tails). As in Papers I and II, the properties of cosmological shock waves are analyzed, the 
energy dissipations to gas thermal energy and CR energy are evaluated, and the vorticity 
distribution is analyzed. We then compare the results for the three simulations to highlight 
the effects of non- gravitational processes on the properties of shocks and their roles on the 
cosmic plasmas in the large scale structure of the universe. 

Simulations are described in §2. The main results of shock identification and properties, 
energy dissipations, and vorticity distribution are described in §3, §4, and §5, respectively. 
Summary and discussion are followed in §6. 



Simulations 



The results reported here are based on the simulations previously presented in lCen fc Ostriker 



( 120061 ). The simulations included radiative processes of heating/cooling, and the two sim- 
ulations with and without galactic superwind (GSW) feedbacks were compared in that pa- 
per. Here an additional adiabatic (non-radiative) simulation with otherwise the same setup 
was performed. Hereafter these three simulations are referred as "Adiabatic", "NO GSW", 
and "GSW" simulations, respectively. Specifically, the WMAPl-normalized ACDM cosmol- 
ogy was employed with the following parameters: Qb = 0.048, Qm = 0.31, Q\ = 0.69, 
h = i^o/(100 km/s/Mpc) = 0.69, erg = 0.89, and n = 0.97. A cubic box of comoving size 
85 h~^Mpc was simulated using 1024^ grid zones for gas and gravity and 512^ particles for 
dark matter. It allows a uniform spatial resolution of Al = 83 /i^^kpc. In Papers I and II, 
an adiabatic simulation in a cubic box of comoving size 100 /i~^Mpc with 1024^ grid zones 
and 512^ particles, employing slightly different cosmological parameters, was used. The 
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simul ations were performed using a PM/Eulerian hydrodynamic cosmology code (iRyu et al 
1993h . 



Detailed descriptions for input physical ingredients such as non-equilibrium ioniza- 
tion/cooling, phot oionization /heat i ng, star formation, an d feedback processes can be found 
in earlier papers (jCen et al\ l2003l : ICen fc Ostrikerl |2006| ). Feedbacks from star formation 
were treated in three forms: ionizing UV photons, GSWs, and metal enrichment. GSWs 
were meant to represent cumulative supernova explosions, and modeled as outflows of sev- 
eral hundred km s~^. The input of GSW energy for a given amount of star formation was 
determined by matching the outflow velocities computed for star- burst gal axies in the sim- 
ulation with those observed in the real world JPettini et a/jl2002h (see also ICen &: Ostriken 
20061 . for details). 



Figure 1 shows the gas mass distribution in the gas density-temperature plane, fmiPgas, T), 
and the gas mass fraction as a function of gas temperature, fm(T), at z = for the three 
simulations. The distributions are quite different, depending primarily on the inclusion of 
radiative cooling and photoionization/heating. GSW feedbacks increase the fraction of the 
WHIM with 10^ < T < 10 ''K, and at the same time affect the distribution of the warm/diffuse 
gas with T < 10^ 



3. Properties of Cosmological Shock Waves 

We start to describe cosmological shocks by briefing the procedure by which the shocks 
were identified in simulation data. The details can be found in Paper I. A zone was tagged as 
a shock zone currently experiencing shock dissipation, whenever the following three criteria 
are met: 1) the gradients of gas temperature and entropy have the same sign, 2) the local 
flow is converging with V ■ v < 0, and 3) |AlogT| > 0.11 corresponding to the temperature 
jump of a shock with M > 1.3. Typically a shock is represented by a jump spread over 
2 — 3 tagged zones. Hence, a shock center was identifled within the tagged zones, where 
V ■ is minimum, and this center was labeled as part of a shock surface. The Mach number 
of the shock center, M, was calculated from the temperature jump across the entire shock 
zones. Finally to avoid confusion from complex flow patterns and shock surface topologies 
associated with very weak shocks, only those portions of shock surfaces with M > 1.5 were 
kept and used for the analysis of shocks properties. 

Figure 2 shows the locations of identifled shocks in a two-dimensional slice at z = in the 
GSW simulation. The locations are color-coded according to shock speed. As shown before 
in Paper I, external accretion shocks encompass nonlinear structures and reveal, in addition 
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to cluster complexes, rich topology of filamentary and sheet-like structures in the large scale 
structure. Inside the nonlinear structures, there exist complex networks of internal shocks 
that form by infall of previously shocked gas to filaments and knots and during subclump 
mergers, as well as by chaotic flow motions. The shock heated gas around clusters extends 
out to ~ 5 h~^Mpc, much further out than the observed X-ray emitting volume. 

In the GSW simulation, with several hundred km s~^ for outflows, the GSW feedbacks 
affected most greatly the gas around groups of galaxies, while the impact on clusters with 
kT > 1 keV was minimal. In Figure 3 we compare shock locations in a region around 
two groups with kT ~ 0.2 — 0.3 keV in the three simulations. It demonstrates that GSW 
feedbacks pushed the hot gas out of groups with typical velocities of ~ 100 km s^^ (green 
points). In fact the prominent green balloons of shock surfaces arou nd groups in Figure 2 



are due to GSW feedbacks (see also Figure 4 of ICen fc Ostrikerl |2006| ) 



In the left panels of Figure 4 we compare the surface area of identified shocks, normalized 
by the volume of the simulation box, per logarithmic Mach number interval, dS{M) / d log M 
(top), and per logarithmic shock speed interval, dSiVs) / dlogVs (bottom), at 2; = in the 
three simulations. Here S and K are given in units of (/i"^Mpc)~^ and km s~^. The 
quantity S provides a measure of shock frequency or the inverse of the mean comoving 
distance between shock surfaces. The distributions of dS{M) / d\ogM for the NO GSW and 
GSW simulations are similar, while that for the Adiabatic simulation is different from the 
other two. This is mainly because the gas temperature outside nonlinear structures is lower 
without photoionization/heating in the Adiabatic simulation. As a result, external accretion 
shocks tend to have higher Mach number due to colder preshock gas. The distribution of 
dSiVs) I d log on the other hand, is similar for all three simulations for > 15 km s~^. For 
Vg < 15 km s~^, however, there are more shocks in the Adiabatic simulation (black points in 
Figure 3). Again this is because in the Adiabatic simulation the gas temperature is colder in 
void regions, and so even shocks with low speeds of Vg < 15 km s~^ were identified in these 
regions. The GSW simulation shows slightly more shocks than the NO GSW simulation 
around Vg ~ 100 km s~^, because GSW feedbacks created balloon-shaped surfaces of shocks 
with typically those speeds (green points in Figure 3). 

For identified shocks, we calculated the incident shock kinetic energy flux, = (l/2)piV^/, 
where pi is the preshock gas density. We then calculated the kinetic energy flux through 
shock surfaces, normalized by the volume of the simulation box, per logarithmic Mach num- 
ber interval, dF^{M) / d\og M , and per logarithmic shock speed interval, dF^iVg) / d log Vg. In 
the right panels of Figure 4, we compare the flux at 2; = in the three simulations. Once 
again, there are noticeable differences in dF^{M)/dlogM between the Adiabatic simulation 
and the other two simulations, which can be interpreted as the result of ignoring photoion- 



- 7- 



ization/heating in the gas outside nonlinear structures in the Adiabatic simulation. GSW 
feedbacks enhance only slightly the shock kinetic energy flux for Vg ~ 100—300 km s~^, as can 
be seen in the plot of dF^{Vs) / dlogVs. Yet, the total amount of the energy flux is expected 
to be quite similar for all three simulations. This implies that the overall energy dissipation 
at cosmological shocks is governed mainly by the gravity of matter, and that the inclusion of 
various non-gravitational processes such as radiative cooling, photoionization/heating, and 
GSW feedbacks have rather minor, local effects. 

We note that a temperature floor of Tfloor = ^cbr was used for the three simulations in 
this work, while Tfloor = lO'' K was set in paper I. It was because in Paper I only an adiabatic 
simulation was considered and the 10^ K temperature floor was enforced to mimic the effect 
of photoionization/heating on the IGM. However we found that when the same temperature 
floor is enforced, the statistics of the current Adiabatic simulation agree excellently with 
those of Paper I. Specifically, the shock frequency and kinetic energy flux, dS{M) / d log M 
and dF^{M) / d\ogM , for weak shocks with 1.5 < M < 3 are a bit higher in the current 
Adiabatic simulation, because of higher spatial resolution. But the total kinetic energy flux 
through shock surfaces, F^{M > 1.5), agrees within a few percent. On the other hand. In 
Paper I we were able to reasonably distinguish external and internal shocks according to the 
preshock temperature, i.e., external shocks if Ti < Tgoor and internal shocks if Ti > Tfloor- We 
no longer made such distinction in this work, since the preshock temperature alone cannot 
tell us whether the preshock gas is inside nonlinear structures or not in the simulations with 
radiative cooling. 



4. Energy dissipation by Cosmological Shock Waves 



The CR injection and acceleration rates at shocks depend in general upon the shock 
Mach number, fleld obliquity angle, and the strength of the Alfven turbulence responsible 
for scattering. At quasi-parallel shocks, in which the mean magnetic fleld is parallel to the 
shock normal direction, small anisotropy in the particle velocity distribution in the local fluid 
frame cau ses some particles in t he high energy tail of the Maxwellian distribution to stream 
upstream ( Giacalone et a/.lll992l ). The streaming motions of the high energy particles against 
the background fluid generate strong MHD Alf ven waves upstream of the shoc k, which in turn 
scatter particles and amphfy magnetic flelds (lBelllll978l : iLucek fc Belli 120001 ). The scattered 
particles can then be a ccelerated further to higher energies via Fermi flrst order process 
( iMalkov fc Druryll200ll ). These processes, i.e., leakage of suprathermal particles into CRs, 
self-excitation of Alfven waves, ampliflcation of magnetic flelds, and further acceleration 
of CRs, are all integral parts of coUisionless shock formation in astrophysical plasmas. It 
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was shown that at strong quasi-parallel shocks, 10~^ — 10"^ of the incoming particles can be 
injected into the CR population, up to 60% of the shock kinetic energy can be transferred into 
CR i ons, and at the same time substantial nonli near feedbacks are exerted to the underlying 
flow JSerezhko et a/.lll995l : E^ang fc Joneslboosi ). 



At perpendicular shocks with weakly perturbed magnetic fields, on the other hand, par- 
ticles gain energy mainly by drifting along the shock surface in the vxB electric fiel d. Such 
drift acceleration can be much more efficient than the accelerat ion at parallel shocks (jJokipii 
19871 : iKang et all 119971 : lOstrowski fc Siemieniec-Oziebld l2002l ). But the particle injection 
into the acceleration process is expected to be inefficient at perpendicu lar shocks, since th e 
transport of particles norrn al to the average field direction is suppressed (lEllison et a/.lll995l ). 
However, iGiacalond (120051 ) showed that the injection problem at perpendicular shocks can be 
alleviated substantially in the presence of fully turbulent fields owing to field line meandering. 

As in Paper I, the gas thermalization and CR acceleration efficiencies are defined as 
6{M) = Fth/F^ and ri{M) = Fcr/F^, respectively, where Fth is the thermal energy fiux 
generated and Fcr is the CR energy fiux accelerated at shocks. We note that for gasdy- 
namical shocks without CRs, the gas thermalization efficiency can be calculated from the 
Rankine-Hugoniot jump condition, as follows: 



5o(M) 



€-th.2 — ^th.l 



V2 



where the subscripts 1 and 2 stand for preshock and postshock regions, respectively. The 
second term inside the brackets subtracts the effect of adiabatic compression occurred at a 
shock too, not just the thermal energy fiux entering the shock, namely, Cth.iVi- 

At CR modified shocks, however, the gas thermalization efficiency can be much smaller 
than Sq{M) for strong shocks with large M, since a significant fraction of the shock kinetic 
energy can be transferred to CRs. The gas thermalization and CR acceleration efficiencies 
were estimated using the results of DSA simulations of quasi-parallel shocks with Bohm 
diffusion coef ficient, self-consisten t treatments of thermal leakage injection, and Alfven wave 
propagation (IKang fc Jones! 120071 ). The simulations were started with purely gasdynamical 
shocks in one-dimensional, plane-parallel geometry, and CR acceleration was followed by 
solving the diffusion-convection equation explicitly with very high resolution. Shocks with 
Vs = 150 — 4500 km s~^ propagating into media of Ti = 10^ — 10^ K were considered. 
After a quick initial adjustment, the postshock states reach time asymptotic values and the 
CR modified shocks evolve in an a pproximately self-sim ilar way with the shock structure 
broadening linearly with time (refer iKang &: JonesI 120071 . for details). Given this self-similar 
nature of CR modified shocks, we calculated time asymptotic values of 6{M) and r]{M) as 
the ratios of increases in the gas thermal and CR energies at shocks to the kinetic energy 
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passed through the shocks at the termination time of the DSA simulations. As in Eq. (1), 
the increase of energies due to adiabatic compression was subtracted. 

Figure 5 shows 6{M) and //(M) estimated from DSA simulations and their fittings for 
the cases with and without a preexisting CR component. The fitting formulae are given 
in Appendix A. Without a preexisting CR component, gas thermalization is more efficient 
than CR acceleration at shocks with M < 5. However, it is likely that weak internal shocks 
propagate through the IGM that contains CRs accelerated previously at earlier shocks. 
In that case, shocks with preexisting CRs need to be considered. Since the presence of 
preexisting CRs is equivalent to a higher inj ection rate, CR accel eration is more efficient in 



that case, especially at shocks with M < 5 (IKang fc Jones II2003I ). In the bottom panel the 
efficiencies for shocks with Pcr/ Pg ~ 0.3 in the preshock region are shown. For comparison, 
(5o(M) for shocks without CRs is also drawn. Both 5{M) and ri{M) increase with Mach 
number, but ri{M) asymptotes to ~ 0.55 while 6{M) to ~ 0.30 for strong shocks with 
M > 30. So about twice more energy goes into CRs, compared to for gas heating, at strong 
shocks. 

The efficiencies for the case without a preexisting CR component in the upper panel 
of Figure 5 can be directly compared with the same quantities presented in Figure 6 of 
Paper I. In Paper I, however, the gas thermalization efficiency was not calculated explicitly 
from DSA simulations, and hence 6o{M) for gasdynamic shocks was used. It represents 
gas thermalization reasonably well for weak shocks with M < 2.5, but overestimates gas 
thermalization for stronger CR modified shocks. Our new estimate for r]{M) is close to that 
in Paper I, but a bit smaller, especially for shocks with M < 30. This is because inclusion of 
Alfven wave drift and dissipation in the shock precu rsor reduces the effect ive velocity change 



experienced by CRs in the new DSA simulations of iKang fc Joned (120071 ). 



A note of caution for ri{M) should be in order. As outlined above, CR injection is 
less efficient and so the CR acceleration efficiency would be lower at perpendicular shocks, 
compared to at quasi-parallel shocks. CR injection and acceleration at oblique shocks are not 
well understood quantitatively. And the magnetic field directions at cosmological shocks are 
not known. Considering these and other uncertainties involved in the adopted DSA model, 
we did not attempt to make further improvements in estimating S{M) and ri{M) at general 
oblique shocks. But we expect that an estimate at realistic shocks with chaotic magnetic 
fields and random shock obliquity angles would give reduced values, rather than increased 
values, for r]{M). So t]{M) given in Figure 5 may be regarded as upper limits. 

By adopting the efficiencies in Figures 5, we calculated the thermal and CR energy fluxes 
dissipated at cosmological shocks, dFth{M) /d log M, dFth{Vs) / dlogVg, dFcR{M)/d\ogM 
and dFcRiVs) /d log Vs, using Fth = F^S{M) and Fcr = F^ri{M), in the same way we 
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calculated dF^{M)/d\ogM and dFff^iVs) / dlogVs in the previous section. We then integrated 
from 2; = 5 to 2; = the shock kinetic energy passed and the thermal and CR energies 
dissipated through shock surfaces as follows: 

dlogX Sth,oJz=5 dhgX 

where the subscript i = (f), th, or CR stands for the kinetic, thermal, or CR energies fluxes, 
the variable X is either M or Vg, and Sth,o is the total gas thermal energy at 2; = inside 
the simulation box normalized by its volume. 

Figure 6 shows the resulting dYi{M)/d\ogM and dYiiVg) / d log Vg and their cumulative 
distributions, Yi{> M) and Yi{> Vg), for the GSW simulation. Weak shocks with M < 4 or 
fast shocks with Vg > 500 km s~^ are responsible most for shock dissipations, as already noted 
in Paper I. While the thermal energy generation peaks at shocks in the range 1.5 ^ M < 3, 
the CR energy peaks in the range 2, 5 < M < 4 if no preexisting CRs are included or in 
the range 1.5 < M < 3 if preexisting CRs of Pcr/ Pg ~ 0.3 in the preshock region are 
included. With our adopted efficiencies, the total CR energy accelerated and the total gas 
thermal energy dissipated at cosmological shocks throughout the history of the universe 
are compared as Ycr{M > 1.5) ~ 0.5Ft/i(M > 1.5), when no preexisting CRs are present. 
With preexisting CRs in the preshock region, the CR acceleration becomes more efficient, 
so Ycr{M > 1.5) ~ 1.7Yth{M > 1.5), i.e., the total CR energy accelerated at cosmological 
shocks is estimated to be 1.7 times the total gas thermal energy dissipated. We note here 
again that these are not meant to be very accurate estimates of the CR energy in the IGM, 
considering the difficulty of modeling shocks as well as the uncertainties in the DSA model 
itself. However, they imply that the IGM and the WHIM, which are bounded by strong 
external shocks with high M and filled with weak internal shocks with low M, could contain 
a dynamically significant CR population. 



5. Vorticity Generation at Cosmological Shock Waves 



Cosmological shocks formed in the large scale structure of the universe are by nature 
curved shocks, accompanying complex, often chaotic flow pa tterns. It is well known that vor- 
ticity , = V X {7, is generated at such curved oblique shocks (lBinneylll974l : iDavies fc Widrow 
2OOOI ). In Paper II, the generation of vorticity behind cosmological shocks and turbulence 
dynamo of magnetic flelds in the IGM were studied in an adiabatic ACDM simulation. In 
this study we analyzed the distribution of vorticity in the three simulations to assess quanti- 
tatively the effects of non-gravitational processes. Here we present the magnitude of vorticity 



- 11 - 



with the vorticity parameter 

r(f, z) ^ t^,.{z)u;{r, z) = :^^y (3) 

where ta.gc{.z) is the age of the universe at redshift z. With teddy = l/f^ interpreted as 
local eddy turnover time, r represents the number of local eddy turnovers in the age of 
the universe. So if r 3> 1, we expect that turbulence has been fully developed after many 
turnovers. 

Figure 7 shows fluid quantities and shock locations in a two-dimensional slice of (21.25 /i^^Mpc)^, 
delineated by a solid box in Figure 2, at 2; = in the GSW simulations. The region contains 
two clusters with kT ~ 1 — 2 keV in the process of merging. Bottom right panel shows that 
vorticity increases sharply at shocks. The postshock gas has a larger amount of vorticity 
than the preshock gas, indicating that most, if not all, of the vorticity in the simulation was 
produced at shocks. 

Figure 8 shows the gas mass distribution in the gas density-vorticity parameter plane, 
fmiPgas, t) , (uppcr panel) and the gas mass fraction per logarithmic r interval, dfm{T) / d log r, 
(bottom panel) for the three simulations. The most noticeable point in the upper panel is 
that vorticity is higher at the highest density regions with p = pgas/ (Pgas) ^ 10^ in the NO 
GSW and GSW simulations than in the Adiabatic simulation. This is due to the additional 
flow motions induced by cooling. Inclusion of GSW feedbacks, on the other hand, does not 
alter signiflcantly the overall distribution in the gas density-vorticity parameter plane. The 
bottom panel indicates that cooling increased the mass fraction with large vorticity r > 10, 
while reduced the mass fraction with 1 < t < 10. GSW feedbacks increased slightly the 
mass fraction with 1 ^ t < 10, which corresponds to the gas in the regions outskirts of 
groups that expand further out due to GSWs {i.e., balloons around groups). But overall we 
conclude that the non-gravitational processes considered in this paper have limited effects 
on vorticity in the large scale structure of the universe. 

We note that the highest density regions in the NO GSW and GSW simulations have 
T ~ 30 on average. As described in details in Paper II, such values of r imply that local 
eddies have turned over many times in the age of the universe, so that the ICM gas there 
has had enough time to develop magnetohydro dynamic (MHD) turbulence. So in those 
regions, magnetic flelds should have grown to have the energy approaching to the turbulent 
energy. On the other hand, the gas with 1 ^ p ^ 10^, mostly in fllamentary and sheet-like 
structures, has 0.1 < r < 10. MHD turbulence should not have been fully developed there 
and turbulence growth of magnetic flelds would be small. Finally in the low density void 
regions with p < 1, vorticity is negligible with r < 0.1 on average, as expected. 
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6. Summary 

We identified cosmological shock waves and studied their roles on cosmic plasmas in 
three cosmological N-body/hydrodynamic simulations for a concordance ACDM universe in 
a cubic box of comoving size 85 h~^Mpc: 1) adiabatic simulation (Adiabatic), 2) simulation 
with radiative cooling and photoionization/heating (NO GSW), and 3) same as the second 
simulation but also with galactic superwind feedbacks (GSW). The statistics and energetics 
of shocks in the adiabatic simulation are in an excellent agreement with those of Paper I 
where an adiabatic simulation with slightly different cosmological parameters in a cubic box 
of comoving size 100 h~^Mpc was analyzed. 

Photoionization/heating raised the gas temperature outside nonlinear structures in the 
NO GSW and GSW simulations. As a result, the number of identified shocks and their 
Mach numbers in the NO GSW and GSW simulations were different from those in the Adi- 
abatic simulation. GSW feedbacks pushed out gas most noticeably around groups, creating 
balloon-shaped surfaces of shocks with speed Vs ~ 100 km s~^ in the GSW simulation. How- 
ever, those have minor effects on shock energetics. The total kinetic energy passed through 
shock surfaces throughout the history of the universe is very similar for all three simula- 
tions. So we conclude that the energetics of cosmological shocks was governed mostly by 
the gravity of matter, and the effects non-gravitational processes, such as radiative cooling, 
photoionization/heating, and GSW feedbacks, were rather minor and local. 

We estimated both the improved gas thermalization efficiency, 5(M), and CR acceler- 
ation efficiency, ri{M), as a function shock Mach number, from nonlinear diffusive shock 
simulations for quasi-parallel shocks that assumed Bohm diffusion for CR protons and in- 
corpora ted self-consistent tr eatments of thermal leakage injection and Alfven wave propa- 



gation (IKang fc JonesI 120071 ). The cases without and with a preexisting CR component of 
Pcr/Pq ~ 0.3 in the preshock region were considered. At strong shocks, both the injection 
and acceleration of CRs are very efficient, and so the presence of a preexisting CR component 
is not important. At shocks with with M > 30, about 55 % of the shock kinetic energy goes 
into CRs, while about 30 % becomes the thermal energy. At weak shocks, on the other hand, 
without a preexisting CR component, the gas thermalization is more efficient than the CR 
acceleration. But the presence of a preexisting CR component is critical at weak shocks, 
since it is equivalent to a higher injection rate and the CR acceleration becomes more effi- 
cient with it. As a result, ri{M) is higher than 6{M) even at shocks with M < 5. However, 
at perpendicular shocks, the CR injection is suppressed, and so the CR acceleration could 
be less efficient than at parallel shocks. Thus our CR shock acceleration efficiency should be 
regarded as an upper limit. 



With the adopted efficiencies, the total CR energy accelerated at cosmological shocks 
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throughout the history of the universe is estimated to be Ycr{M > 1.5) ~ 0.5 Yth{M > 1.5), 
i.e., 1/2 of the total gas thermal energy dissipated, when no preexisting CRs arc present. 
With a preexisting CR component of Pcn/Pg ~ 0.3 in the prcshock region, Ycr{M > 1.5) ~ 
1.7 Yth{M > 1.5), i.e., the total CR energy accelerated is estimate to be 1.7 times the total 
gas thermal energy dissipated. Although these are not meant to be very accurate estimates of 
the CR energy in the ICM, they imply that the ICM could contain a dynamically significant 
CR population. 

We also examined the distribution of vorticity inside the simulation box, which should 
have been generated mostly at curved cosmological shocks. In the ICM, the eddy turn-over 
time, teddy = l/*-^) is about 1/30 of the age of the universe, i.e., r = tage/ teddy ~ 30. In 
filamentary and sheet- like structures, r ~ 0.1 — 10, while r < 0.1 in void regions. Radiative 
cooling increased the fraction of gas mass with large vorticity r > 10, while reduced the 
mass fraction with 1 < r < 10. GSW feedbacks increased shghtly the mass fraction with 
1 < r < 10. Although the effects of these non-gravitation effects are not negligible, the 
overall distribution of vorticity arc similar for the three simulations. So we conclude that the 
non-gravitational processes considered in this paper do not affect significantly the vorticity 
in the large scale structure of the universe. 
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NNG05GK10G and NSF grant AST-0507521. The work of HK and DR was also supported in 
part by Korea Foundation for International Cooperation of Science & Technology (KICOS) 
through the Cavendish-KAIST Research Cooperation Center. 



The gas thermalization efficiency, 5(M), and the CR acceleration efficiency, f]{M), for 
the case without a preexisting CR component (in upper panel of Figure 5) are fitted as 
follows: 
for M < 2 



A. Fitting Formulae for S{M) and r]{M) 



S{M) = 0.92 So 
r]{M) = 1.96 X 10-^(M2 - 1) 



(Al) 
(A2) 



for M > 2 



4 



(M-1) 



n 




(A3) 



n=0 
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ao = -4.25, ai = 6.42, aa = -1.34, 03 = 1.26, 04 = 0.275 

n=0 

bo = 5.46, 61 = -9.78, 62 = 4.17, 63 = -0.334, 64 = 0.570 



(A4) 
(A5) 
(A6) 



The efficiencies for tfie case witfi a preexisting CR component (in bottom panel of Figure 
5) are fitted as follows: 



for M < 1.5 



for M > 1.5 



S{M) = 0.90 So 
r]{M) = 1.025 So 



(A7) 
(A8) 

(A9) 



n=0 



ao = -0.287, ai = 0.837, 03 = -0.0467, 03 = 0.713, 04 = 0.289 (AlO) 

^ . 1/1// — I I'" 

(All) 



n=0 



bo = 0.240, 61 = -1.56, 62 = 2.80, 63 = 0.512, 64 = 0.557 



(A12) 



Here So{M) is the gas thermalization efficiency at shocks without CRs, which was cal- 
culated from the Rankine-Hugoniot jump condition, (black solid line in Figure 5): 



So{M) = 



7(7 - l)M2i? 



- (7 - 1) 
, (7 + 1) 

7 + 1 



pi 7 - 1 + 2/M2 



(A13) 
(A14) 
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Fig. 1. — Top panels: Gas mass distribution in the gas density-temperature plane at z = 
for the Adiabatic, NO GSW, and GSW simulations. Bottom panel: Gas mass fraction as a 
function of gas temperature at z = for the three simulations. 




Fig. 2. — Two-dimensional slice of (85 h^^Mpc)^ showing shock locations at z = in 
the GSW simulation, which are color-coded according to shock speed as follows: black for 
< 15 km s~^, blue for 15 < < 65 km s~^, green for 65 < < 250 km s~^, red for 
250 < Vs < 1000 km s~^, and magenta for Vg > 1000 km s^^. A blown-up image of the box 
(dashed line) in the upper right corner is shown in Figure 3, while a blown-up image of the 
box (solid line) around two merging clusters is shown in Figure 7. 
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Fig. 3. — Two-dimensional slice of (21.25 /i^^Mpc)^ showing shock locations at z = in 
the Adiabatic, NO GSW and GSW simulations. The locations are color-coded according to 
shock speed. Two groups in the GSW simulation have kT ~ 0.2 — 0.3 keV. 
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Fig. 4. — Left panels: Inverse of the mean comoving distance between shock surfaces as a 
function of Mach number M (top) and shock speed Vg (bottom) ai z = for the Adiabatic 
(solid line), NO GSW (dashed line), and GSW (dotted line) simulations. Right panels: 
Kinetic energy flux per comoving volume passing through shock surfaces in units of 10^° ergs 
s~^ {h~^Mpc)~^ as a function of M (top) and Vg (bottom). Note that the bottom two panels 
have different ranges of abscissa. 




Fig. 5. — Gas thermalization efficiency, S{M), and CR acceleration efficiency, ri{M), as 
a function of Mach number. Red and blue dots are the values estimated from numerical 
simulations based on a DSA model and red and blue lines are the fits. The top panel shows 
the case without preexisting CRs, while the bottom panel shows the case with preexisting 
CRs at a level of Pcr/Pq ~ 0.3 in the preshock region. Black solid line is for the gas 
thermalization efficiency for shocks without CRs. 
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Fig. 6. — Left panels: Shock kinetic energy passed, dY^ (dotted line), thermal energy dis- 
sipated, dYth (dashed line), and CR energy dissipated, dYcR (solid line), through surfaces 
of cosmological shocks with Mach number between logM and logM + d{\ogM) (top) and 
through surfaces of cosmological shocks with shock speed between log Vg and log Vs+d{\og Vs) 
(bottom), integrated from z = 5 to z = 0. Red and magenta lines are the CR energy for the 
cases without and with preexisting CRs, respectively. Blue and green lines are the thermal 
energy for the cases without and with preexisting CRs, respectively. The thermal energy 
expected to be dissipated at cosmological shocks without CRs (long dashed cyan line) is also 
plotted for comparison. Right panels: Cumulative energy distributions, 1^(> M) (top) and 
Yi{> Vs) (bottom), for Mach number greater than M and for shock speed greater than Vs. 
The energies are normalized by the gas thermal energy at z = inside the simulation box 
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Fig. 7. — Two-dimensional slice of (21.25 /i^^Mpc)^ around two merging clusters with kT ~ 
1 — 2 keV at 2; = in the GSW simulation. Distributions of gas density (top left), temperature 
(top right), shock locations (bottom left), and vorticity (bottom right) are shown. In the 
gas density, temperature, and vorticity distributions, back, blue and red contours represent 
regions of low, middle, and high values, respectively. 
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Fig. 8. — To'p panels: Gas mass distribution in the gas mass density- vorticity parameter 
plane at 2; = for the Adiabatic, NO GSW, and GSW simulations. The vorticity parameter 
is defined as r = utage{z), where u = \V x v\ and tage{z) is the age of the universe at redshift 
z. Bottom panel: Gas mass fraction as a function of vorticity parameter at z = for the 
three simulations. 



